#include "Molmodel.h"
#include "Simbody.h"
#include <fstream>
#include <ostream>
#include <cstdlib>
#include <cmath>
#include <exception>
#include <string.h>

using namespace SimTK;
using namespace std;
const char* PdbFileName="pdboutput.pdb";
const string Output_folder="E:\\Sorted Files\\User Files\\Documents\\Research summer 2012\\Sim\\output\\";

class OutputReporter : public PeriodicEventReporter {
	
public:
    OutputReporter(const MultibodySystem& system, const SimbodyMatterSubsystem& matter, const NoseHooverThermostat& Therm, ostream& PosFile, 
		ostream& VelFile, ostream& AccFile, ostream& K_Energy, ostream& P_Energy, ostream& Temp, Real reportInterval) 
    :    PeriodicEventReporter(reportInterval), system(system), matter(matter), Therm(Therm), PosFile(PosFile), VelFile(VelFile), AccFile(AccFile), K_Energy(K_Energy), P_Energy(P_Energy), Temp(Temp) {
		
	}
	/*
    ~OutputReporter()
	{
		PosFile.close();
		VelFile.close();
		AccFile.close();
		Energy.close();
		Temp.close();
	}
	*/
	void handleEvent(const State& state) const{
       //system.realize(state, Stage::Position);
		
		PosFile<<state.getTime()<<"\t";
		VelFile<<state.getTime()<<"\t";
		AccFile<<state.getTime()<<"\t";
        for(int i=0; i<state.getNQ(); i++)
		{
			PosFile<<state.getQ()[i]<<"\t";
			VelFile<<state.getQDot()[i]<<"\t";
			AccFile<<state.getQDotDot()[i]<<"\t";

		}
		K_Energy<<state.getTime()<<"\t"<<matter.calcKineticEnergy(state);
		P_Energy<<state.getTime()<<"\t"<<system.calcPotentialEnergy(state);
		Temp<<state.getTime()<<"\t"<<Therm.getCurrentTemperature(state);

		PosFile<<endl;
		VelFile<<endl;
		AccFile<<endl;
		K_Energy<<endl;
		P_Energy<<endl;
		Temp<<endl;
			
    }


private:
    const MultibodySystem&  system;
	const SimbodyMatterSubsystem& matter;
	const NoseHooverThermostat& Therm;
	ostream& PosFile;
	ostream& VelFile;
	ostream& AccFile;
	ostream& K_Energy;
	ostream& P_Energy;
	ostream& Temp;

};

int main() {
try{
	//molecule-specialized simbody System
	CompoundSystem system;
	SimbodyMatterSubsystem matter(system);
	GeneralForceSubsystem forces(system);
	DecorationSubsystem decorations(system);

	//molecular force field
	DuMMForceFieldSubsystem forceField(system);
	forceField.loadAmber99Parameters();


	int equilib_sample_size = 10;
	int steptime=30;

	ofstream pdbFile;
	ofstream PosFile;
	ofstream VelFile;
	ofstream AccFile;
	ofstream K_Energy;
	ofstream P_Energy;
	ofstream Temp; 
	
	pdbFile.open(Output_folder+PdbFileName);
	PosFile.open(Output_folder+"Positions.txt", ios::app);
	VelFile.open(Output_folder+"Velocities.txt", ios::app);
	AccFile.open(Output_folder+"Accelerations.txt", ios::app);
	K_Energy.open(Output_folder+"K_Energies.txt", ios::app);
	P_Energy.open(Output_folder+"P_Energy.txt", ios::app);
	Temp.open(Output_folder+"Temps.txt", ios::app);

	//forceField.setCoulombGlobalScaleFactor(0);
	//forceField.setGbsaGlobalScaleFactor(0);

	RNA rna("AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA");
	rna.assignBiotypes();
	

	NoseHooverThermostat thermo(forces, matter, 300,0.07,6);

	system.addEventReporter(
		new Visualizer::Reporter(system, 0.050));

	
	system.addEventReporter(
		new PeriodicPdbWriter(system, pdbFile, 0.2500));

	system.addEventReporter(
		new OutputReporter(system, matter, thermo, PosFile, VelFile, AccFile, K_Energy, P_Energy, Temp, 0.25));

	rna.setCompoundBondMobility(BondMobility::Free);
	/*
	for(int i=0;i<rna.getNumResidues();i++){
		rna.setResidueBondMobility(ResidueInfo::Index(i),BondMobility::Rigid);
	}
	*/
	system.adoptCompound(rna);
	system.modelCompounds();

	State state = system.realizeTopology();

	LocalEnergyMinimizer::minimizeEnergy(system, state, 15.0);
	cout<<"num residues"<<rna.getNumResidues();
	
	VerletIntegrator integ(system, .001);
	TimeStepper ts(system, integ);
	ts.initialize(state);
	ts.stepTo(steptime);

	PosFile.close();
	VelFile.close();
	AccFile.close();
	K_Energy.close();
	P_Energy.close();
	Temp.close(); 
	pdbFile.close();
	std::cout << "\nWrote pdb file "
              << Pathname::getCurrentWorkingDirectory()
			  << PdbFileName << "\n";
    return 0;
} 
catch(const std::exception& e) {
    std::cerr << "ERROR: " << e.what() << std::endl;
    return 1;
}
catch(...) {
    std::cerr << "ERROR: An unknown exception was raised"
              << std::endl;
    return 1;
}
}